home *** CD-ROM | disk | FTP | other *** search
/ Aminet 33 / Aminet 33 - October 1999.iso / Aminet / misc / math / TCalcStats2c.lha / TCalcStats2c / AREXX / anova_oneway.rexx next >
Encoding:
OS/2 REXX Batch file  |  1998-09-22  |  11.1 KB  |  588 lines
  1. /* ANOVA - One Way */
  2.  
  3. options results
  4. if ~show('P','TCALC') then do
  5.     address command 'run turbocalc:turbocalc'
  6.     address command 'waitforport TCALC'
  7.     loadflag=1
  8. end
  9. address 'TCALC'
  10. 'DEFPUBSCREEN()'
  11. /* Add-in Rexx Libraries needed for some routines */
  12. signal on syntax
  13. if ~show('l','rexxmathlib.library') then
  14.    call addlib('rexxmathlib.library',0,-30)  /* add to library list */
  15. if ~show('l','rexxreqtools.library') then
  16.    call addlib('rexxreqtools.library',0,-30)
  17. if ~show('l','rexxsupport.library') then
  18.    call addlib('rexxsupport.library',0,-30)
  19. signal off syntax
  20.  
  21. /* Start Main Routine */
  22. if loadflag=1 then 'Load()'
  23. 'ActivateWindow()'
  24. range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  25. colon=pos(":",range)
  26. if colon=0 then do
  27.     'Message "Please select a range before executing this script"'
  28.     'DEFPUBSCREEN("Workbench")'
  29.     exit
  30. end
  31.  
  32. /* Find cell references and cell, column numbers */
  33. start_cell=substr(range,1,colon-1)
  34. end_cell=substr(range,colon+1)
  35. start_row=cellrow(start_cell)
  36. end_row=cellrow(end_cell)
  37. start_col=cellcol(start_cell)
  38. end_col=cellcol(end_cell)
  39. NRows=end_row-start_row+1
  40. NCols=end_col-start_col+1
  41.  
  42. /* Get cell reference for output range */
  43. out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request")/*,,'rt_pubscrname="TCALC"')*/
  44. if out_cell="" then exit
  45. if length(out_cell)<2 | datatype(left(out_cell,1),'n')=1 then do
  46.     'Message "Invalid cell reference"'
  47.     'DEFPUBSCREEN("Workbench")'
  48.     exit
  49. end
  50.  
  51. /* Suppress Screen Redraw to Speed Things Up */
  52. 'Refresh 0'
  53.  
  54. /* Open a small output window on tcalc screen*/
  55. fo=0
  56. CR='0a'x
  57. DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
  58. if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
  59.      call writeln(6Info, DisplayMsg)
  60.     fo=1
  61. end
  62. else do
  63.     'Message "TCALC Screen not available for Progress messages"'
  64. end
  65. CALL DELAY(150)
  66.  
  67. /* Get cell references for top cell in each column */
  68. 'SelectCell' start_cell
  69. do col=start_col to end_col
  70.     'GetCursorPos'
  71.     top_cell.col=result
  72.     'Column' 1
  73. end
  74.  
  75. /* Get labels for later use on output */
  76. 'SelectCell' start_cell
  77. 'GetValue'
  78. testlabel=result
  79. testlabel=strip(testlabel)
  80. if datatype(testlabel,'n')=1 then do
  81.     labelflag=0
  82.     do x=1 to NCols
  83.     title.x="Column "||x
  84.     end
  85. end
  86. else do
  87.     labelflag=1
  88.     NRows=NRows-1
  89.     do x=1 to NCols
  90.     'GetValue'
  91.     str=result
  92.     title.x=translate(strip(str),"_"," ")
  93.     'Column 1'
  94.     end
  95. end
  96. if fo then call writech(6Info,"Progress...10 ")
  97.  
  98. /* Get data from cell range */
  99. col=start_col
  100. lav=0
  101. tot=0
  102. count.=0
  103. total.=0
  104. do x=1 to NCols
  105.     'SelectCell' top_cell.col
  106.     if labelflag=1 then 'CursorDown 1'
  107.     do y=1 to NRows
  108.         'GetValue'
  109.         valtest=result
  110.         if datatype(valtest)='NUM' then do
  111.             'GetValue'
  112.             val=result
  113.             data.x.y=val
  114.             tot=tot+val
  115.             total.x=tot
  116.             count.x=1+count.x
  117.         end
  118.         'CursorDown 1'
  119.     end
  120.     col=col+1
  121.     tot=0
  122.     lav=0
  123.     val=0
  124. end
  125.  
  126. if fo then call writech(6Info,"20 ")
  127.  
  128. /* Calculate Means */
  129. mean.=0
  130. do x=1 to NCols
  131.     mean.x=total.x/count.x
  132. end
  133. if fo then call writech(6Info,"30 ")
  134.  
  135. /* Calculate Standard deviation and Variance */
  136. dat=0
  137. meenx=0
  138. sum.=0
  139. sum3.=0
  140. sum4.=0
  141. sd.=0
  142. var.=0
  143. m2.=0
  144. m3.=0
  145. m4.=0
  146. do x=1 to NCols
  147.     sum.x=0
  148.     sum3.x=0
  149.     sum4.x=0
  150.     meenx=mean.x
  151.     do y =1 to count.x
  152.     dat=data.x.y
  153.     sum.x=(dat-meenx)**2+(sum.x)
  154.     sum3.x=(dat-meenx)**3+(sum3.x)
  155.     sum4.x=(dat-meenx)**4+(sum4.x)
  156.     end
  157.     N=(count.x)-1
  158.     var.x=(sum.x)/N
  159.     sd.x=sqrt(var.x)
  160.     m2.x=(sum.x)/(count.x)
  161.     m3.x=(sum3.x)/(count.x)
  162.     m4.x=(sum4.x)/(count.x)
  163. end
  164.  
  165. if fo then call writech(6Info,"40 ")
  166.  
  167. /* Calculate Totals */
  168. TN.=0
  169. sumTN=0
  170. GrandTotal=0
  171. NumObs=0
  172. TNG=0
  173. Do x=1 to NCols
  174.     NumObs=NumObs+count.x
  175.     GrandTotal=GrandTotal+total.x
  176.     TN.x=((total.x)**2)/count.x
  177.     SumTN=SumTN+TN.x
  178. end
  179. TNG=(GrandTotal**2)/NumObs
  180. if fo then call writech(6Info,"50 ")
  181.  
  182. /* Calculate Sum of Squares */
  183. SS.=0
  184. SSTotal=0
  185. SSBN=0
  186. SSWN=0
  187. SST=0
  188. do x=1 to NCols
  189.     do y=1 to count.x
  190.         SS.x=((data.x.y)**2)+SS.x
  191.     end
  192.     SSTotal=SSTotal+SS.x
  193. end
  194. if fo then call writech(6Info,"60 ")
  195. SSBN=SumTN-TNG /* Sum of Squares Between */
  196. SSWN=SSTotal-SumTN /* Sum of Squares Within */
  197. SST=SSTotal-TNG /* Total Sum of Squares */
  198.  
  199. /* Calculate Variances */
  200. DFW=NumObs-NCols /* Degrees of Freedom Within */
  201. DFB=NCols-1 /* Degrees of Freedom Between */
  202. DFT=NumObs-1 /* Degress of Freedom for Total */
  203. VarB=SSBN/DFB /* Variance Estimate for Between */
  204. VarW=SSWN/DFW /* Variance Estimate for Within */
  205. VarT=VarB/VarW /* F Ratio for Total */
  206. if fo then call writech(6Info,"70 ")
  207.  
  208. /* Calculate Probability */
  209.  
  210. AC=.0001
  211. EC=.005
  212. EC2=EC+EC
  213. P=PROBABILITY(VarT,DFB,DFW)
  214. P=1-P
  215. if fo then call writech(6Info,"80 ")
  216. /* P=.95 OR .99 */
  217. FCRIT1=FCRITICAL(.95,DFB,DFW)
  218. if fo then call writech(6Info,"90 ")
  219. FCRIT2=FCRITICAL(.99,DFB,DFW)
  220. if fo then do
  221.     call writeln(6Info,"100 ")
  222.     call writeln(6Info,"Writing output to window...")
  223. end
  224. /* Output */
  225. 'SelectCell' out_cell
  226. 'ColumnWidth 15'
  227. 'Put "ANOVA - One Way"'
  228. 'CursorDown 2'
  229. 'Put' "Group"
  230. 'Column 1'
  231. 'ColumnWidth 10'
  232. 'Alignment 2'
  233. 'Put' "Count"
  234. 'Column 1'
  235. 'ColumnWidth 10'
  236. 'Alignment 2'
  237. 'Put' "Total"
  238. 'Column 1'
  239. 'ColumnWidth 12'
  240. 'Alignment 2'
  241. 'Put' "Mean"
  242. 'Column 1'
  243. 'ColumnWidth 10'
  244. 'Alignment 2'
  245. 'Put' "Variance"
  246. 'SelectCell' out_cell
  247. 'CursorDown 3'
  248. do x=1 to NCols
  249.     'GetCursorPos'
  250.     first_cell.x=result
  251.     title=""""||title.x||""""
  252.     'Put' title
  253.     'CursorDown 1'
  254. end
  255.  
  256. Do x=1 to NCols
  257.     'SelectCell' first_cell.x
  258.     'Column 1'
  259.     'Put' count.x
  260.     'Column 1'
  261.     'Put' total.x
  262.     'Column 1'
  263.     'Put' format(mean.x,,4)    
  264.     'Column 1'
  265.     'Put' format(var.x,,4)
  266. end
  267.  
  268. 'SelectCell' first_cell.NCols
  269. 'CursorDown 1'
  270. 'CursorDown 1'
  271. 'Put "ANOVA"'
  272. 'CursorDown 1'
  273. 'Put "Source of Variation"'
  274. 'GetCursorPos'
  275. Go_cell=result
  276. 'Column 1'
  277. 'Alignment 2'
  278. 'Put "SS"'
  279. 'Column 1'
  280. 'Alignment 2'
  281. 'Put "d.f. (v):"'
  282. 'Column 1'
  283. 'Put "Variance Est.:"'
  284. 'Column 1'
  285. 'Alignment 2'
  286. 'Put "F-Ratio:"'
  287. 'SelectCell' Go_cell
  288. 'CursorDown 2'
  289. 'Put "Between Groups:"'
  290. 'CursorDown 1'
  291. 'Put "Within Groups:"'
  292. 'CursorDown 1'
  293. 'Put "Total:"'
  294. 'CursorDown 2'
  295. 'Put "P(F<=f):"'
  296. 'CursorDown 1'
  297. 'Put "F-Critical One-Tail(95%):"'
  298. 'CursorDown 1'
  299. 'Put "F-Critical One-Tail(99%):"'
  300. 'SelectCell' Go_cell
  301. 'CursorDown 2'
  302. 'Column 1'
  303. 'Put' format(SSBN,,4)
  304. 'CursorDown 1'
  305. 'Put' format(SSWN,,4)
  306. 'CursorDown 1'
  307. 'Put' format(SST,,4)
  308. 'CursorDown 2'
  309. 'Put' format(P,,6)
  310. 'CursorDown 1'
  311. 'Put' format(FCRIT1,,4)
  312. 'CursorDown 1'
  313. 'Put' format(FCRIT2,,4)
  314. 'CursorUp 6'
  315. 'Column 1'
  316. 'Put' DFB
  317. 'CursorDown 1'
  318. 'Put' DFW
  319. 'CursorDown 1'
  320. 'Put' DFT
  321. 'CursorUp 2'
  322. 'Column 1'
  323. 'Put' format(VarB,,4)
  324. 'CursorDown 1'
  325. 'Put' format(VarW,,4)
  326. 'CursorUp 1'
  327. 'Column 1'
  328. 'Put' format(VarT,,4)
  329. 'Refresh 1'
  330. 'Refresh 2'
  331. /*'Message' "Finished"*/
  332. /*indicate the main script is finished*/
  333. DisplayMsg="Cleaning up ...."||CR||"Exiting"
  334. result=writeln(6Info, DisplayMsg)
  335. if result~=0 then do
  336.     /*Wait 3 seconds*/
  337.     CALL DELAY(150)
  338.     /* close window*/
  339.     result=close(6Info)
  340. end
  341. 'DEFPUBSCREEN("Workbench")'
  342. exit
  343.  
  344. /* Procedures */
  345.  
  346. cellrow: procedure
  347. do
  348.     parse arg cell
  349.     do charpos=2 to length(cell)
  350.     if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
  351.     end
  352.     return 0
  353. end
  354. Return
  355.  
  356. cellcol: procedure
  357. do
  358.     parse arg cell
  359.     labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
  360.     cell=upper(cell)
  361.     len=length(cell)
  362.     val=0
  363. do charpos=1 to len
  364.     if datatype(substr(cell,charpos,1),n) then
  365.     do cell=reverse(substr(cell,1,charpos-1))
  366.     do x=1 to length(cell)
  367.     val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
  368.     end
  369.     return val
  370.     end
  371.     end
  372.     return 0
  373. end
  374. Return
  375.  
  376. syntax:
  377.      if arg(1)='FAIL' then do
  378.         'Message "Library is unavailable."'
  379.         'DEFPUBSCREEN("Workbench")'
  380.         exit
  381.         end
  382.     'Message "Unknown Syntax Error"'
  383.     'DEFPUBSCREEN("Workbench")'
  384.     exit
  385.  
  386. Format:  procedure
  387.  
  388.      arg number, before, after
  389.      CallLine = SIGL
  390.      if ~datatype(CallLine, 'N') then CallLine = '??'
  391.  
  392.      /* Make sure we have a number as first (required) argument    */
  393.      if ~datatype(number, 'N') then do
  394.         if number = '' then
  395.            fc = 17     /* Wrong number of arguments           */
  396.         else
  397.            fc = 47     /* Arithmetic conversion error             */
  398.         signal FormatSyntaxError
  399.      end
  400.      num = number + 0
  401.      if before = '' & after = '' then
  402.         return num
  403.      else do
  404.         parse var num integer '.' fraction
  405.         if before = '' then before = length(integer)
  406.         if after = '' then after = length(fraction)
  407.         if ~datatype(before, N) | ~datatype(after, N) then
  408.            do fc = 18
  409.            signal FormatSyntaxError
  410.        end
  411.         if before < length(integer) then do
  412.            fc = 18
  413.            signal FormatSyntaxError
  414.         end
  415.         if after ~= length(fraction) then do
  416.            fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
  417.         if integer<1&integer>-1 then integer=integer
  418.            else integer = integer + (fraction % 1)
  419.            fraction = substr(fraction, 3)
  420.         end
  421.         if fraction >= 0 then
  422.            return right(integer, before)'.'fraction
  423.         else
  424.            return right(integer, before)
  425.      end
  426.  
  427.  FormatSyntaxError:
  428.         if show('F', STDERR) then
  429.            call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
  430.         else
  431.            mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
  432.         'Message' mess
  433.         parse source Func .
  434.         if Func = 'FUNCTION' then do
  435.         'DEFPUBSCREEN("Workbench")'
  436.            exit "Err"
  437.         end
  438.         else do
  439.         'DEFPUBSCREEN("Workbench")'
  440.            exit 10
  441.         end
  442.  
  443. FCRITICAL: PROCEDURE EXPOSE AC EC EC2
  444.  
  445.     PARSE ARG PO,K1,K2
  446.     /* fIND NORMAL Z CORRESPONDING TO P */
  447.     P1=PO
  448.     Z=NORMALPP(PO)
  449.     IF K2<3 THEN DO
  450.         W=K2**.75
  451.         W=Z/W
  452.         Z=Z*(1.1581-W*(.2296+W*(.0042+.0027*W)))
  453.     END
  454.     /* FIND INITIAL APPROX. TO F */
  455.     A1=2/(9*K1)
  456.     A2=2/(9*K2)
  457.     W=Z*Z
  458.     W1=1+A2*(A2-W-2)
  459.     W2=A1+A2-A1*A2-1
  460.     W3=1+A1*(A1-W-2)
  461.     SN=SIGN(W2*W2-W1*W3)
  462.     IF SN=-1 THEN W3=-1*SQRT(ABS(W2*W2-W1*W3))
  463.     ELSE W3=SQRT(W2*W2-W1*W3)
  464.     FO=(W3-W2)/W1
  465.     FO=FO**3
  466.     /* MODIFIED NEWTON ITERATION TO IMPROVE VALUE OF F */
  467.     DO FOREVER
  468.         FCRIT=FO
  469.         PO=PROBABILITY(FCRIT,K1,K2)
  470.         F1=PO-P1
  471.         FCRIT=FO+EC
  472.         PO=PROBABILITY(FCRIT,K1,K2)
  473.         F2=PO
  474.         FCRIT=FO-EC
  475.         PO=PROBABILITY(FCRIT,K1,K2)
  476.         F2=F2-PO
  477.         F2=F2/EC2
  478.         F1=FO-F1/F2
  479.         IF ABS(F1-FO)>AC THEN FO=F1
  480.         ELSE LEAVE
  481.     END
  482.     FCRIT=F1
  483. RETURN FCRIT
  484.  
  485. NORMALPP: PROCEDURE
  486.  
  487.     ARG P0
  488.     A1=2.515517
  489.     A2=.802853
  490.     A3=.010328
  491.     A4=1.432788
  492.     A5=.189269
  493.     A6=.001308
  494.     Q=.5-ABS(P0-.5)
  495.     SN=SIGN(-2*LN(Q))
  496.     IF SN=-1 THEN W=-1*SQRT(ABS(-2*LN(Q))
  497.     ELSE W=SQRT(-2*LN(Q))
  498.     W1=A1+W*(A2+A3*W)
  499.     W2=1+W*(A4+W*(A5+A6*W))
  500.     ZZ=W-W1/W2
  501.     SELECT
  502.         WHEN (P0-.5)<0 THEN TT=-1
  503.         WHEN (P0-.5)=0 THEN TT=0
  504.         otherwise TT=1
  505.     END
  506.     ZZ=ZZ*TT
  507. RETURN ZZ
  508.  
  509. PROBABILITY: PROCEDURE EXPOSE AC EC EC2
  510.  
  511.     PARSE ARG F,K1,K2
  512.     IF K1=1 THEN AC=.00001
  513.     H2=K2/2
  514.     H1=K1/2
  515.     H3=H1+H2
  516.     L1=0
  517.     XX=H1
  518.     LG=LOGGAMMA(XX)
  519.     L1=L1+LG
  520.     XX=H2
  521.     LG=LOGGAMMA(XX)
  522.     L1=L1+LG
  523.     XX=H3
  524.     LG=LOGGAMMA(XX)
  525.     L1=L1-LG
  526.     W1=K2/(K2+K1*F)
  527.     XX=0
  528.     IF H2<(H3*W1) THEN DO
  529.         W2=W1
  530.         W1=1-W1
  531.         W3=H2
  532.         H2=H1
  533.         H1=W3
  534.     END
  535.     ELSE DO
  536.         W2=1-W1
  537.         XX=1
  538.     END
  539.     T1=1
  540.     W3=1
  541.     P=1
  542.     I=INT(H1+W2*H3)
  543.     W4=W1/W2
  544.     /*LABELB:*/
  545.     T2=H1-W3
  546.     IF I=0 THEN W4=W1
  547.     /*LABELC:*/
  548.     DO FOREVER
  549.         T1=T1*T2*W4/(H2+W3)
  550.         P=P+T1
  551.         T2=ABS(T1)
  552.         IF (T2<=AC) & (T2<=AC*P) THEN LEAVE /*SIGNAL 'LABELD'*/
  553.         W3=W3+1
  554.         I=I-1
  555.         IF I>=0 THEN DO /*SIGNAL 'LABELB'*/
  556.             T2=H1-W3
  557.             IF I=0 THEN W4=W1
  558.             ITERATE
  559.         END
  560.         T2=H3
  561.         H3=H3+1
  562.     /*SIGNAL 'LABELC'*/
  563.     END
  564.     /*LABELD:*/
  565.     W1=EXP(H2*LN(W1)+(H1-1)*LN(W2)-L1)
  566.     P=P*W1/H2
  567.     IF XX=0 THEN PO=P
  568.     ELSE PO=1-P
  569.     AC=.001
  570. RETURN PO
  571.  
  572. LOGGAMMA: PROCEDURE
  573.  
  574.     ARG XA
  575.     C1=76.18009173
  576.     C2=-86.50532033
  577.     C3=24.01409822
  578.     C4=-1.231739516
  579.     C5=.001208580
  580.     C6=-.000005364
  581.     C7=2.506628275
  582.     X1=XA-1
  583.     W=X1+5.5
  584.     W=(X1+.5)*LN(W)-W
  585.     S=1+C1/(X1+1)+C2/(X1+2)+C3/(X1+3)
  586.     S=S+C4/(X1+4)+C5/(X1+5)+C6/(X1+6)
  587.     L=W+LN(C7*S)
  588. RETURN L